SoSe2021

Folienübersicht

Neue Komponenten.

..die es bei mehr als 1 erklärenden Variablen zu berücksichtigen gilt:

  1. Zusätzliche Datenexploration
  2. Modellauswahl
  3. Zusätzliche Modellvalidierung
  4. Anpassung der Modellgüte
  5. Interpretation des numerischen R Output und Dummy-Variablen
  6. Visualisierung der Ergebnisse

Alle X kontinuierlich: Multiple Regression

Multiple Regression

Die Prinzipien der einfachen Regressionsanalyse können auf zwei oder mehr unabhängige Variablen erweitert werden:


Faustregel: nicht mehr als 10 unabhängige Variablen (X)

Voraussetzungen für multiple lineare Regressionsmodelle

  1. Unabhängigkeit
  2. Homogenität / homogene Varianz
  3. Normalität / Normalverteilung
  4. Linearität

Neu

  1. Anzahl der Kovariaten < Anzahl der Beobachtungen (p < n)
  2. Keine lineare Abhängigkeit zwischen den unabhängigen Variablen (keine exakte Multikollinearität)

Problem der Multikollinearität

  • Schätzungen von \(\beta\)s und Signifikanztests unzuverlässig.
    • Kleine Änderungen in den Daten führen zu extremen Änderungen in den Regressionskoeffizienten.
    • Schätzungen zwischen verschiedenen Stichproben der gleichen Grundgesamtheit können sich stark unterscheiden (→ hohe Standardfehler der geschätzten Parameter).

Prüfung

  • Paarweise Streudiagramme (mit pairs()).
  • Berechnung der Korrelationsmatrix (mit cor()).
  • Zusätzliche Bestimmung der Toleranz (TOL) und des VIF-Werts (= ‘Variance Inflation Factor’) für jede Variable.

Variance Inflation Factor

Kurz: VIF

\[VIF = \frac{1}{TOL} = \frac{1}{1-R^2}\]

  • \(R^2\) hier aus Modell welches alle X Variablen enthält außer die, für die der VIF-Wert berechnet wird.
  • Faustregel: VIF > 3 wird als problematisch eingestuft (manche sagen VIF > 10, hängt auch von der Anzahl der unabhängigen Variablen ab)
  • VIF in R: car::vif()
    • Wenn X Variablen VIF > 3 haben, Variable mit höchstem VIF entfernen (oder nach biologischem Wissen entscheiden),
    • VIFs für verbleibende Variablen neu berechnen,
    • und fortfahren, bis alle Variablen VIF < 3 haben.

Variance Inflation Factor

Simulation

Angenommen y ist hauptsächlich abhängig von x1 und nun soll der Steigungskoeffizient \(\beta_1\) geschätzt werden → wenn andere X Variablen, die von x1 abhängen, mit modelliert werden, variiert die Schätzung für \(\beta_1\) bei jeder Wiederholung sehr stark:

Funktionen

# Funktionen
calc_coefs <- function(x1, x2, x3){
  y <- x1 + rnorm(length(x1), sd = .3)
  # beta für x1 aus 3 Modellen extrahieren
  c(coef(lm(y ~ x1))[2],  
    coef(lm(y ~ x1 + x2))[2],  
    coef(lm(y ~ x1 + x2 + x3))[2]) 
}

run_sim <- function(x1, x2, x3, 
  n_sim = 1000) {
  betas <- sapply(1:n_sim, 
    function(i) calc_coefs(x1, x2, x3))
  # Berechnung der Varianz
  round(apply(betas, 1, var), 5)
}

Vergleich

# Zufällige X Variablen erzeugen
n <- 100
x1 <- rnorm(n)
x2a <- rnorm(n)
x3a <- rnorm(n)
# x2 und x3 in Abhängigkeit zu x1:
x2b <- x1/sqrt(2) + rnorm(n) /sqrt(2)
x3b <- x1*0.95 + rnorm(n)*sqrt(1-0.95^2)

run_sim(x1, x2a, x3a)
     x1      x1      x1 
0.00099 0.00099 0.00099 
run_sim(x1, x2b, x3b)
     x1      x1      x1 
0.00097 0.00229 0.01117 

Schätzen und Testen der Koeffizienten

  • Auch hier erfolgt die Schätzung über die OLS Methode.
  • Y-Achsenabschnitt (Konstante) \(\alpha\):
    • Wert von Y wenn alle X Variablen Null betragen → ökologisch weniger relevant
  • Partielle Regressionskoeffizienten (Steigungen) \(\beta_p\):
    • Einfluss auf Y wenn \(X_p\) sich um eine Einheit erhöht, wobei alle anderen unabhängigen Variablen oder Prädiktoren konstant gehalten werden.
    • Separate t-Tests für jeden partiellen Regressionskoeffizienten im Modell: \(H_0: \beta_p = 0\)

Anpassung des Bestimmtheitsmaß

(Multiple) \(R^2\)

\[R^2 = \frac{SS_{Regression}}{SS_{Gesamt}} = 1-\frac{SS_{Residuen}}{SS_{Gesamt}}= 1-\frac{\sum(y_i-\hat{y}_i)^2}{\sum(y_i-\bar{y})^2}\]

  • Problem bei multiplen Regression:
    • Mit zunehmender Zahl an unabhängigen Variablen steigt die Gesamtstreuung im Nenner, auch wenn diese nicht zur Erklärung beitragen.
    • \(\sum(y_i-\bar{y})^2\rightarrow \infty~~~\Rightarrow R^2 \rightarrow 1\)

Adjusted \(R^2\)

\[R_{adj}^2 =1-(1-R^2)\left(\frac{(N-1)}{(N-p-1)}\right)\]

  • → d.h. die Freiheitsgrade werden mit berücksichtigt (p = Anzahl der unabhängigen Variablen)

Wie vergleicht und wählt man das beste Modell aus?

  • Basierend auf der besten Anpassungsgüte (‘goodness-of-fit’): Wie gut kann ein Modell die Beobachtungen vorhersagen?
  • Wie beurteilen wir das?
    • Visuell
    • Verwendung eines Kriteriums, dass die Diskrepanz zwischen beobachteten Werten und den zu erwartenden Werten des Modells zusammenfasst, z.B.
      • erklärte Varianz (\(R^{2}\))
      • ‘Akaike’s Information Criterion’ (AIC) und Modifikationen
      • ‘Bayesian Information Criterion’ (BIC)
  • Ein gutes Buch für diese Thematik wurde von Anderson & Burnham im Jahr 2002 geschrieben: Model Selection and Multimodel Inference - A Practical Information-Theoretic Approach

‘Akaike’s Information Criterion’

AIC()

  • Der AIC ist eine relative Messgröße, keine absolute wie \(R^2\).
    • Die Werte können nur zwischen Modellen verglichen werden, die an die gleichen Beobachtungen und die gleiche y Variable angepasst wurden.
    • →Du kannst keine Modelle vergleichen, die für verschiedene Teilmengen angepasst wurden oder wenn y in Teilen transformiert wurde!
  • Je kleiner der AIC, desto besser.
  • Faustregel: Wenn der AIC Unterschied zwischen 2 Modellen < 2 ist, dann wähle das einfachere Modell.
  • Bei großen Stichproben empfiehlt sich das Bayesian Information Criterion (BIC), da hier der Strafterm ab 8 Variablen ein größeres Gewicht erhält als beim AIC.

Beispiel: airquality

Haben Wind und Sonneneinstrahlung einen Einfluss auf die in New York gemessenen Ozon-Werte?

str(airquality)
'data.frame':   153 obs. of  6 variables:
 $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
 $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
 $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
 $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
 $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
 $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...
# counts of NA
purrr::map_int(airquality, ~sum(is.na(.)))
  Ozone Solar.R    Wind    Temp   Month     Day 
     37       7       0       0       0       0 
air <- airquality %>% tidyr::drop_na(Ozone) %>% tidyr::drop_na(Solar.R)

Beispiel: airquality

Datenexploration

hist(air$Ozone)

air$OzoneS <- sqrt(air$Ozone)
hist(air$OzoneS)

pairs(air[ ,c(7,2,3)])

car::vif(lm(OzoneS ~ Wind + Solar.R, data = air))
    Wind  Solar.R 
1.016442 1.016442 

Beispiel: airquality

Mögliche Modelle

# Nullmodell:
mod0 <- lm(OzoneS ~ 1, data = air)

mod1 <- lm(OzoneS ~ Wind, data = air)
mod2 <- lm(OzoneS ~ Solar.R, data = air)
mod3 <- lm(OzoneS ~ Wind + Solar.R, data = air)

# Maximalmodell:
mod4 <- lm(OzoneS ~ Wind + Solar.R + Wind:Solar.R, data = air)

Beispiel: airquality

Vergleich des \(R^2\)

summary(mod1)$r.squared
summary(mod2)$r.squared
summary(mod3)$adj.r.squared # Da 2 X, hier nun das angepasste R^2 auswählen
summary(mod4)$adj.r.squared
[1] 0.3703109
[1] 0.162264
[1] 0.4682859
[1] 0.4846256
  • Wenn beide erklärenden Variablen im Modell sind, wird deutlich mehr Variabilität erklärt.
  • Das Maximalmodell mit der Interaktion bringt aber kaum ein Gewinn an erklärter Variabilität.

Beispiel: airquality

Vergleich des AIC

AIC(mod0, mod1, mod2, mod3, mod4)
     df      AIC
mod0  2 516.0493
mod1  3 466.7085
mod2  3 498.3965
mod3  4 447.8994
mod4  5 445.4023
  • Das komplexeste Modell hat den niedrigsten AIC.
  • Allerdings beträgt die Differenz zum weniger komplexen Modell nur 2.5 → Wahl fällt auf mod3.

‘All subsets selection’

  • Die eben vorgestellte Methode der Modellauswahl nennt sich auch ‘all subsets selection’.
  • Hier werden alle möglichen Teilmengen von Modellen verglichen und das optimale Modell anschl. ausgewählt.
  • Nachteile:
    • Kann zu einer großen Anzahl von zu vergleichenden Modellen führen.
    • Viele Iterationen und zeitaufwendig, wenn p > 5.
  • Die Anzahl der Modelle kann jedoch reduziert werden, indem a priori einige wenige zu vergleichende Modelle ausgewählt werden.
  • Einige bevorzugen diesen Ansatz, da er einige sehr gute Modelle hervorbringen kann.

Beispiel: airquality

Modellvalidierung 1 - mod3

par(mfrow = c(2,2))
plot(mod3)

Beispiel: airquality

Modellvalidierung 2 - mod3

Zusätzliche Plots: Residuen vs. jede X Variable

air$residuals <- resid(mod3)
par(mfrow = c(1,2), mar = c(4,4,1,1))
plot(x = air$Wind, y = air$residuals); abline(h = 0, lty = 2)
plot(x = air$Solar.R, y = air$residuals); abline(h = 0, lty = 2)

Beispiel: airquality

Ergebnis

summary(mod3)
Call:
lm(formula = OzoneS ~ Wind + Solar.R, data = air)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.9188 -1.3407 -0.1783  1.4120  3.9330 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.241963   0.647547  12.728  < 2e-16 ***
Wind        -0.388544   0.048079  -8.081 9.96e-13 ***
Solar.R      0.008855   0.001877   4.719 7.14e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.779 on 108 degrees of freedom
Multiple R-squared:  0.478, Adjusted R-squared:  0.4683 
F-statistic: 49.44 on 2 and 108 DF,  p-value: 5.706e-16

Beispiel: airquality

Visualisierung - Code

Verwende nicht die geom_smooth() Funktion!!

library(modelr)
p_w <- air %>% data_grid(
  Wind = seq_range(Wind, 20), # 20 Werte im gleichen Abstand zwischen min und max
  Solar.R = mean(Solar.R) # auf Null setzen oder Konstante wie MW wählen
) %>% 
  add_predictions(mod3) %>%
  ggplot(aes(x = Wind)) + ylab("OzoneS") +
  geom_line(aes(y = pred), color = "red") +
  geom_point(data = air, mapping = aes(y = OzoneS)) 
p_s <- air %>% data_grid(
  Wind = mean(Wind), # diese nun konstant halten
  Solar.R = seq_range(Solar.R, 20)
) %>% 
  add_predictions(mod3) %>%
  ggplot(aes(x = Solar.R)) + ylab("OzoneS") +
  geom_line(aes(y = pred), color = "blue") +
  geom_point(data = air, mapping = aes(y = OzoneS)) 
gridExtra::grid.arrange(p_w, p_s, ncol = 2)

Beispiel: airquality

Visualisierung - Plot

Fragen?